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Abstract 

Efficient eiectronic structure methods can be built around efficient tensor representations of the wavefunc- 
tion. Here we describe a general view of tensor factorization for the compact representation of electronic 
wavefunctions. We use these ideas to construct low-complexity representations of the doubles amplitudes in 
local second order M0ller-Plesset perturbation theory. We introduce two approximations - the direct orbital 
specific virtual approximation and the full orbital specific virtual approximation. In these approximations, 
each occupied orbital is associated with a small set of correlating virtual orbitals. Conceptually, the repre- 
sentation lies between the projected atomic orbital representation in Pulay-Saeb0 local correlation theories 
and pair natural orbital correlation theories. We have tested the orbital specific virtual approximations on 
a variety of systems and properties including total energies, reaction energies, and potential energy curves. 
Compared to the Pulay-Saeb0 ansatz, we find that these approximations exhibit favourable accuracy and 
computational times, while yielding smooth potential energy curves. 
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I. INTRODUCTION 



In electron correlation there are two problems of complexity. The first relates to the information 
(storage) required to represent the wavefunction and the second to the complexity of manipulating 
the wavefunction to calculate observables. Consider, for example, the doubles amplitudes tfj 
(where ij denote occupied orbitals, ab virtual orbitals) common to the second-order M0ller-Plesset 
perturbation (MP2), coupled cluster doubles (CCD), and coupled electron pair approximations 
(CEPA). The storage scales like iV 4 while the cost of obtaining the energy scales like N 5 for MP2 
and iV 6 for CCD and CEPA, where N is a measure of the size of the system. 

In several limits, we expect these formal complexities to be too high. For example, if there is a 
large number of atomic orbitals on a single center, there is redundancy in the product occupied- 
virtual pair basis. Also, in large molecules the cost to obtain the energy should be linear in 
the size of the molecule. In both these situations, the mismatch between formal complexity and 
our expected complexity suggests that the amplitudes and amplitude equations have some special 
structure. For example, in Pulay-Saeb0 local correlation theories based on the projected atomic 
orbital (PAO) ansatz [IH2] the sparsity structure of the amplitudes is built in through distance- 
based truncations [4], and using this framework it has been possible to achieve linear scaling of 
storage and computational cost with system size [SHS] . Naturally, these truncations are not orbitally 
invariant and require a representation of the orbitals in which the amplitude matrix is maximally 
sparse. An important task in devising algorithms with reduced complexities is to find optimal 
transformations of the amplitudes to representations that are approximable with low complexity. 

The doubles amplitude tfj is a tensorial quantity. Consequently, the problem of finding a low 
complexity representation can be viewed as one of tensor representation or tensor factorization. 
While there has been much recent work in constructing low complexity representations for the dou- 
bles amplitudes and two-electron integrals, both through using more optimal orbitals (for example, 
along the lines of pair natural orbitals (PNO) [1UH16] . optimized virtual orbital spaces |17l 118] . 
frozen natural orbitals [19-22], and others |23| I24j) as well as matrix factorizations of the integrals 
and amplitudes (such as Cholesky decompositions (CD) and density fitting (DF) (or resolutions 
of the identity, RI) ) . these approaches have not yet explored the full generality of the 

tensorial structure. Exploring the possibility of a more general class of tensor factorizations is the 
concern of the current work. 

Here we first start with a general classification of the different kinds of tensor factorizations. 
Next we introduce a specific factorization which we refer to as an orbital specific virtual approx- 
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imation. This representation has a simple and intuitive interpretation that bridges earlier work 
on optimal virtuals and work on pair natural orbitals. We explore the orbital specific virtual ap- 
proximation in the context of local second order M0ller-Plesset perturbation theory. We present 
benchmark applications on a variety of large molecules, clusters, and reactions. We find that the 
ansatz is very favourable both in its formal properties such as potential energy curve smoothness, 
and weak computational dependence on the size of the underlying basis, as well as regarding its 
absolute costs in terms of storage and timings when compared to an existing efficient implementa- 
tion of the local Pulay-Saeb0 correlation ansatz [6j [36]. Finally, we finish with a discussion of the 
future prospects of such an approach. 

II. THEORY AND ALGORITHM 

A. Classification of tensor factorizations 

We are concerned primarily with the doubles amplitude tensor tfj. We illustrate it pictorially as 
a connected four-point object (see objects on the left in Figs. [I] and [2]). A closely related quantity, 
particularly in second-order M0ller-Plesset theory, is the two-electron integral vfj. In canonical 
closed-shell MP2 theory, the two are related by 



To construct a low- complexity representation of the two-electron integrals or amplitudes, we 
approximate the high-dimensional amplitude or integral tensors by lower-dimensional components. 
These component tensors may share the same "physical" indices i,j,a,b as the target tensor, 
but may also carry additional "auxiliary" indices \,[A,i/,p,. . .. Since the auxiliary indices do not 
appear in the target tensor, they must be traced over in some way, and both the distribution 
of the physical and auxiliary indices amongst the component tensors, as well as the pattern of 
contractions, defines the particular tensor representation. Since these contractions are usually 
non- linear, it is often useful to visualize the contractions pictorially rather than algebraically. 

To illustrate this, consider first the density- fitting and Cholesky decomposition approximations. 
In all these approximations, the two-electron integrals are viewed as a matrix factorization as 



A 

where A is the auxiliary index. Pictorially, we view the above as separating the ia, jb electron-hole 
degrees of freedom, which must then be reconnected via an auxiliary index (see Fig. [T^,) . Naturally, 
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(1) 




(2) 



the efficiency of the factorization relies on the rank of the decomposition (the number of terms in 
the sum) being low. As another example, consider the types of correlation ansatz (such as the 
Pulay-Saeb0 local correlation ansatz) which use a non-canonical virtual orbital basis, for example 
the projected atomic orbital (PAO) virtuals. Non-canonical virtuals are related to canonical 
virtuals (j) a by a transformation 

^ = E*>« ( 3 ) 

and consequently, the canonical and non-canonical doubles amplitudes are related by 

*? = £*W£ ( 4 ) 

which defines the amplitude approximation. Pictorially, this approximation is illustrated in Fig.[2^i. 
For appropriate (such as defined by the PAO virtuals), representation Q allows one to favourably 
exploit locality. For example, in the Pulay-Saeb0-Werner-Schiitz approaches [IrO EH2] , a sparsity 
structure on t^- is imposed by requiring = when i is far apart from j, and for other ij, the 
sum over fx and v is restricted to defined domains [ij] that are in the spatial vicinity of ij. 

Thus, the essence of low-complexity tensor approximation is captured by the types of indices 
on the components and their connectivity, as illustrated in their pictorial representation. We can 
consider generalizations of the above approximations in a variety of ways. For example, we can 
consider approximations with additional auxiliary indices. One example is (cf. Fig. [TJd) 

+ab _ \ "» J.A1A2J.C2 + b J.A4A1 /r\ 

l ij ~ 2^ l i t A 2 A 3 t A 3 A4 t j W 
A1...A4 

= tr[t i t a t 6 t i ] (6) 

where the amplitude tfj is reconstructed as a trace of a matrix product. This approximation recalls 
the matrix product factorization in the density matrix renormalization group (DMRG) [37} 138] . 
although the physical content here is quite distinct, since the DMRG is carried out in the occupation 
number space rather than in the excitation space. Another way to construct new approximations 
is to introduce components with repeated physical indices. A particularly simple example is 

t-ij — tijtiatjb (7) 

where there are no auxiliary indices at all. This recalls the correlator product state approximation 
(also known as an entangled plaquette state [391 HO]), although once again the tensor is expressed 
in an excitation rather than occupation number picture. Naturally, we can consider many other 
combinations of auxiliary indices and physical indices, and the appropriateness of the particular 
choice depends on the problem at hand. 



B. Orbital specific virtual approximation 



We now consider a simple tensor factorization of the doubles amplitudes that we will study in 
this work. We first define the direct orbital specific virtual (dOSV) approximation to the amplitudes 
as (cf. Fig. |J;) 

This has a simple physical interpretation: the component t^ (and similarly t^ J b ) defines a set of 
virtual orbitals for each occupied orbital, and t^ Vj represents amplitudes in this orbital specific 
basis. Note that the subscript i in jUj is somewhat redundant, but we retain it to emphasize that 
fii labels an orbital specific virtual associated with occupied orbital i. By choosing a good set 
of components tf*, either by direct optimization or otherwise (see later) we may define suitable 
adaptations of the virtual basis for each occupied orbital. This is quite natural in a local correlation 
theory, as the optimal orbital specific virtuals for a localized occupied orbital must be located in 
close spatial proximity; however, even when the occupied orbitals are delocalized, we can still 
expect this factorization to be beneficial, as a given occupied orbital does not correlate equally 
with all parts of the virtual space. 

In the direct orbital specific virtual approximation, occupied orbital % excites only to its virtual 
set Hi (i — > fj,i), and occupied orbital j only to its orbital set Vj (j — > Uj), the "exchange" excitations 
% — > — > fii being excluded. It was shown, however, in the context of Pulay-Saeb0 local theory, 
that the inclusion of exchange excitations can lead to greatly improved results |41H44j . While 
formally the exchange excitations can be included by increasing the size of sets \ii and Vj, we can 
also include them explicitly in the structure of the ansatz, which leads to the full orbital specific 
virtual (OSV) approximation 

.ab _ \ ~* f+IH+HiViM , ±V>i + Vi v j + Vj Mj jPjVi Ui , .Mj + P-i v i + v i\ iq\ 

Hj — \ ia L ij L jb L ia L ij L jb "•" L ia L ij L jb "> L ia L ij L jb J \ v > 

fiu 

This can be written in matrix form 

tf = Y( iff* & ) hj H > > b (10) 

^ v V f ij f ij ) V t jb J 

To gain further understanding, we briefly discuss the connection to other approximations where 

non-canonical virtual spaces are used. In theories which use a global set of non-canonical virtuals, 

such as the projected atomic orbital virtual space in the Pulay-Saeb0 ansatz, the amplitude tensor 
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is expressible as in Eq. (j4J), where and t b v parametrize the transformation to a new virtual basis. 
In pair natural orbital theories pioneered some time ago |10H14j and which have been recently 
revisited by Neese and coworkers \Vb\ HE] , a correlating virtual space is defined for each occupied 
pair ij, and the amplitudes are factorized as (cf. Fig. [2)3) 

j.ab _ \ ~* iWj JH v j 

Hj - 2^ '<>'</ l ijb- ' 1 1 > 

Compared to the Pulay-Saebo ansatz where a global set of virtuals is used, the orbital specific 
virtual approximation is able to adapt the virtual space, which leads to a more compact repre- 
sentation of the amplitudes. On the other hand, the orbital specific virtual approximation adapts 
each virtual space to a single occupied orbital rather than a pair. This avoids some of the complex- 
ities inherent to the pair natural orbital ansatz where the definitions of the virtual spaces involve 
four-index components t^ 3 a of similar formal complexity to the amplitudes themselves, and which 
lead to complicated overlap and Fock matrices in the pair natural orbital virtual blocks. Conse- 
quently, we see that formally the orbital specific virtual approximation interpolates between the 
Pulay-Saebo form and the pair natural orbital approximation. We now turn towards its practical 
implementation in second-order perturbation theory (MP2). 



III. IMPLEMENTATION 

A. MP2 wavefunction and singular value orbital specific virtuals 

The central task of MP2 theory is to determine the first-order wavefunction l^^ 1 ^), 

ijab 

where in the above i, j, ■ ■ ■ and a, b, ■ ■ ■ refer respectively to the occupied and virtual spatial orbitals. 
Explicitly in spin-orbital notation 

l*#>= E l^v) (13) 

a,(T'G{Q,^} 

For closed-shell molecules the spin-free orbital notation avoids the explicit use of spin coordinates 
and is very convenient. 



Inserting the orbital specific virtual approximation (cf. Eq. (18])) into Eq. (12), we parametrize 



the first-order wavefunction in terms of factorized amplitudes t^ Uj and corresponding determinants 
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3 °). For the direct orbital specific virtual approximation of MP2 (dOSVMP2), 



i* (1) ) = ^ECi*D ( 14 ) 

where 

while for the full orbital specific virtual approximation of MP2 (OSVMP2) we have 

\jl (>■;" n" ■ c ^f> ■ <t # r> ■ <P ( 16 ) 

ijliu 



In Eq. (14) we have introduced an orbital specific virtual (OSV) defined through an orbital 
transformation from the virtual \a) using tf*, 

lw> = £Ck>. (i7) 



Naturally, we would like the OSVs to be well adapted to each occupied orbital in Eq. (17). One 
quick and economical scheme to determine tf^ is to perform a singular value decomposition (SVD) 
of the MP2 diagonal amplitudes tf^ for each occupied orbital i, 

*sf = £cw&- (is) 

where s Mi is the singular value. In the canonical basis, is directly calculated from, 

tt = Vf " • (19) 

When the localized occupied orbitals are used, ej is the diagonal element of the Fock matrix in the 
local orbital basis. e a and e& are the diagonal elements of the virtual block of Fock matrix. 

The SVD provides a natural setting to truncate the OSV space, as the singular vectors with small 
singular values s M should contribute little to the final amplitudes. Consequently it is reasonable 
to include only those vectors with the largest eigenvalues, keeping either a fixed number of OSVs 
per occupied orbital, a fixed percentage of OSVs, or by using a numerical threshold on s M . We 
have used the first two truncation schemes in this work. After truncation, the complete virtual 
space is parametrized by an incomplete set of orbital specific virtuals. This of course introduces 
errors relative to canonical MP2 theory. However, as numerically shown in the next section, the 
resulting correlation energies exhibit only minor deviations (e.g. < 0.01%) from canonical values, 
while achieving very substantial gains in computational efficiency. 



The OSVs defined above are not always orthogonal. The overlap matrix between the OSVs 
and \uj) is, 

S^ = (^) = Y,C^i (20) 

a 

Through the SVD, the OSVs belonging to the same occupied orbital are orthonormal but the 
OSVs from different occupied orbitals are not. 

We also note that the SVD does not necessarily yield the most optimal orbital specific virtual 
orbitals. The direct optimization of relaxes the OSVs and may help achieve a more compact 
description of the correlation effects. This is under investigation and will be presented elsewhere. 



B. Residual equations 

We derive the exact MP2 residual equations in the orbital specific virtual basis starting from 
the Hylleraas functional, 

h = (*W|F - #(°) IttW) + 2(¥^\Y\^), (21) 

where F is the Fock operator and V is the two-electron fluctuation potential, respectively. is 
the sum of occupied Hartree-Fock eigenvalues. By parametrizing Vl/W in the direct orbital specific 



virtual approximation (dOSVMP2, cf. Eq. (14)) and making the first derivative of h with respect 
to Uj vanish, we arrive at the following formal residual, 

rt = =($ri F - £(o) i* (1) )+( $ ri v i* (o) )=°- ( 22 ) 

The expansion of R^- Uj has the following explicit form for a particular pair 



R (i,j) ~ 2K (i,j) ~ C (hj) + 2 



T (hj) F Cjj) + F (i,i) T (iJ) ~ Yl F kj T (i,k)S(k,j) ~ Yl F ik S (i,k)T(k,j) 



T (j,i) F (i,j) ~ X/ jPifcT (i^) S 



(k,i) 



S(jj). (23) 



Here T(£j) is the matrix form of amplitudes {t^ 3 }, fixing i,j. Fik and F^j are elements of the 
occupied block of the Fock matrix. ^-Uj) and Cuj\ denote the matrices storing two-electron 
Coulomb and exchange integrals ((pi^^il^j^uj) and |0j</> w ), respectively. Suj-\ and Fuj) are, 

respectively, the overlap and Fock matrices with elements (^J^,-) and ((j) i _ Li \F\(j) Uj ) for a pair 
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None of the above quantities gives rise to real bottlenecks for memory or disk storage if the orbital 
specific virtual space is truncated. 



In the full orbital specific virtual approximation (OSVMP2, cf. Eq. (16)), we minimize the 



Hylleraas functional with respect to the amplitudes tfj u % 3 , *, t^ 3 giving rise to analogues 



of ( 22 ) . The explicit residuals are then 



%J) - K (M)+X^ S (w fc ) T (^) [ 6 kj F (ik,ij) ~ FkjS(ik,ij)] + [ S ikF(ij,kj) ~ F ik$(ij,kj)] T fcj S (fcj,ij)} • 
k 

(24) 

where R(ij) is now a matrix of dimension dim(yUj) + dim(i^), with elements of types 
R^- v% , R^ 3 , R^- Ul , Rfj Uj , and similarly for f£(ij)- S(ij,ik) ^ s the overlap matrix between the bras 
{{fjti\, an d the kets \vk)}, and Fuj^fy is the analogously defined Fock matrix. 

C. Projective residual equations for dOSVMP2 (dOSVMP2-P) 

Compared to the standard local MP2 residual equations, those for the orbital specific virtuals 
appear more complicated, especially for the direct (dOSVMP2) ansatz. Formally, we can obtain 
the MP2 amplitudes not only through variational minimization of the Hylleraas functional, but 
also by projection with an appropriate set of bra states. In the case of dOSVMP2, this leads to a 
simpler set of residual equations which define a different set of amplitudes than those arising from 



Eq. (|23|). We write 

jjw>i = ^p\Y - #(o)|,p(i)) + ($g^|V|#(°)) = 0. (25) 

The bra states are chosen to be bi-orthonormal to the ket states in the spatial orbital basis following 
Ref. [HHE], 

(ZiPKi 11 ) = H \ x WjiS^ + SjkSuS^) (26) 

with a normalization prefactor. The bi-orthonormal OSV bra state {^^ Vj \ is defined from the 
canonical bi-orthonormal bra state ($^-|, 

(*iTl = E<*?lOX (27) 

ab 

with [H] 

w ' = 1 r 

1 



(*"l = o( 2 (*tfl + (*i 



(*#|*#) = \ (S lk 5 jl 5 ac 5 bd + 5 jk S u S bc S ad ) (28) 
V 1 + Oijd 
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Note that if we truncate the OSV space, then the space spanned by the bras is not the 

same as the space spanned by {(^- V] |} , due to the presence of the exchange-like excitation ($^| 
in the definition of ($^|. This is what gives rise to the difference between the projective and 
(standard) Hylleraas based residual equations for the dOSVMP2 approximation. 



Thus expanding Eq. (25) yields the projective residual equation R(ij), 



%j) - K (i,j) + T (i,j) F UJ) + f (v) T (m) - Yl F kj T (i,k)S(k,j) ~ Yl F ikS(i,k)T( k ,j) (29) 



which has clearly a simpler form than the standard residual in Eq. ( 23 ) . We denote the direct orbital 
specific virtual approximation defined by the amplitudes from (29), the projective dOSVMP2-P 
approximation. In a complete virtual space, the solution of either dOSVMP2 or dOSVMP2-P 
yields the exact canonical MP2 energy and the first-order wavefunction. In an incomplete virtual 
space (e.g. if not all the orbital specific virtuals are used), however, dOSVMP2-P generally does 
not result in a correlation energy that is variationally bounded above the canonical MP2 value with 
respect to the size of incomplete virtual space. The numerical comparison between dOSVMP2 and 
dOSVMP2-P will be given in the next section. 

D. Preconditioning 

As the Fock matrix is diagonal in the canonical basis, the canonical MP2 amplitudes are directly 
calculated as, 

t$ = 11 , (30) 

where e's are the diagonal elements of the canonical Fock matrix. However, in the orbital specific 
virtual space, the Fock matrix contains significant off-diagonal elements and the residual equations 
must be solved iteratively. Here the use of a preconditioner is essential. To this end we define 
pseudo- virtual energies. The pseudo- virtual energies can be obtained by diagonalizing the Fock 
matrix in the space of orbital specific virtuals for each diagonal occupied pair (£,£), similarly to as 
done in Ref. [5J 

= S(j i j)X(j ji )E(j i j), (31) 

where X(^) is the transformation matrix that diagonalizes both F^) and Su{\. E^n is diagonal 
and contains the pseudo-virtual energies. When preconditioning the residual R(ij) corresponding 
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to orbital pair we then consider the residual as a matrix with the first index corresponding 

to the virtuals associated with i, and the second with the virtuals associated with j. Consequently, 
the transformation is performed on the virtual indices of the original residual R(ij) with Xt ... and 
X(j j) , respectively. 

R (i,3) =X \i,i) R (iJ) X UJ)- ( 32 ) 

The amplitude update matrix At(jj) in the transformed basis is then, 



Al^' = - , (33) 

c l T c 1 c «, c £< 



where ' is the element of the residual matrix R(ij). The pseudo- virtual energies of and e v . 
are the elements of ^ and Eyj) , respectively. Finally the update in the original orbital specific 
basis is given by the back-transformation 

At (M) = X (i,i) Al (i,j) X \j,j)- ( 34 ) 

The same algorithm can be used with both dOSVMP2 and OSVMP2, the latter requiring the 
obvious generalizations of indices to U,j)- 

Note that eigenvectors Xuj\ may be linearly dependent. (The same issue arises e.g. in the 
Pulay-Saeb0 local correlation theory [5]). The redundant vectors are eliminated by a canonical 
orthogonalization of discarding eigenvalues below a given threshold. Here we use a threshold 
of 1(T 6 . 

With the above preconditioning the amplitudes and MP2 energies converge very quickly. For 
example, typically, the change in correlation energy falls below 10 -6 a.u. within 8 ~ 10 iterations 
for dOSVMP2-P and OSVMP2 and within 16 ~ 20 iterations for dOSVMP2. 



E. Computational cost and screening 



The most expensive term in our current implementation of the residual equations ( 23 ) and ( 29 ) 
is the contraction ^ fc FkjTufySfaj) which has a formal scaling of 3 V' i , i.e. iV 6 with system size. 
This is to be contrasted with the canonical MP2 scaling of iV 5 . The local Pulay-Saeb0 ansatz has a 
similar contraction, but because there is a single set of underlying virtual orbitals, the contraction 
can be implemented in 2 V 3 + 5 V 2 operations, which is still iV 5 scaling. In large molecules, V 
is independent of molecular size, and so the scaling of local OSVMP2 is O 3 where the O 3 term has 
a larger prefactor than in the Pulay-Saeb0 theory. 
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We can lower the computational cost by introducing some screening approximations to the iV 6 
contraction. As we have discussed previously, the most important OSVs that can be correlated 
with each occupied orbital are located closely around the occupied orbital itself. By exploiting 
orbital locality, the overlap matrix element S^ kVj decays exponentially with the separation of k 
and j. Based on this we can simply ignore entire classes of unimportant iV 6 contractions in the 
residuals without losing much accuracy, and consequently reduce the formal scaling of solving the 
residual equations in a large system to order O 2 . The following ratio for a given pair (k,j) is 
computed in order to define an appropriate screening threshold, 



IkVk 



= (35) 



With a screening threshold of Tg we then neglect any overlap matrix $(kj\ belonging to a pair (k,j) 
if tjy < Tg. From the definition t| ■ ranges between and 1: thus when Tg = the overlap matrices 
belonging to all pairs of (k,j) are taken into account; when Tg = 1 only the diagonal contributions 



with k = j are included. Eq. (35) allows us to avoid a less desirable spatial truncation criterion. 

Currently, however, the main cost in our OSVMP2 implementation is the integral transforma- 
tion, since the full local occupied space is employed throughout transforming complete two-electron 
integrals. For example, the first-quarter integral transformation (a/3|7<5) — > (if3\j5), which scales 
as N 5 , limits the size of molecules that we can treat efficiently. As has been previously demon- 
strated by Werner and coworkers (5j [6], however, a linear scaling transformation algorithm can be 
achieved by discarding spatially distant occupied orbital pairs and exploiting integral prescreening 
techniques. DF/RI techniques in local approximations [34, [35] can further decrease the cost of 
integral transformations by 1-2 orders of magnitude for large molecules. These techniques will be 
incorporated into our algorithm in future work. 



IV. COMPARING THE dOSVMP2 AND OSVMP2 APPROXIMATIONS 

We have introduced two related factorizations of MP2 theory, the direct orbital specific virtual 
approximation (dOSVMP2) and the full orbital specific virtual approximation (OSVMP2). In 
addition, in dOSVMP2 we can define the amplitudes through two different residual equations: one 
obtained via the Hylleraas functional (dOSVMP2), and one obtained by projection (dOSVMP2- 
P). We now assess the numerical behaviour of these different schemes for correlation energies and 
reaction energies. 
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A. Correlation energies 

The number of OSVs per occupied orbital (the number of \ii or Vi in e.g. Eq. Q) needed 
to recover a given accuracy in the correlation energy relative to the canonical MP2 energy is 
reported for polyglycine oligopeptides, water clusters and polyene chains. As seen in Table |IJ 
a small number of OSVs recovers most of the correlation energy (e.g., > 99.5%), the precise 
number depending on the electronic structure of the molecule. The full OSVMP2 (which includes 
the exchange excitations) requires far fewer OSVs to reach the same accuracy than the direct 
dOSVMP2, typically less than half. Nonetheless both approximations are very compact. The 
number of OSVs to reach a given accuracy also becomes independent of the total molecular size 
very rapidly. For OSVMP2 the correlation energy is saturated at accuracies of 99.5%, 99.9% and 
99.99%, respectively, with 20, 29 and 47 orbitals for [gly] n and with 13, 18 and 29 orbitals for 
(H20)„. Note that this saturation behaviour is expected when the system size becomes much 
larger than its correlation length. 

Compared to [gly] n and (H\20) n , the polyene molecule exhibits significant electronic derealiza- 
tion, and this leads to longer correlation lengths and more extended orbitals. As a result, it is more 
difficult to converge the correlation energy towards the canonical limit than in other molecules. 
For example, with 40 OSVs the OSVMP2 error increases from 0.46% to 0.66% (cf. Table |Tl} as 
the length of polyene chain increases from CgHs to C14H16. Note that this decrease in accuracy is 
physical and not a failure of extensivity of the theory: as the HOMO-LUMO gap of the polyenes 
decreases with increasing chain length, the correlation length increases. 

Regarding the different residual equations for the direct dOSVMP2 factorization, both yield 
very similar results. Using the Hylleraas residual (dOSVMP2), we need a few more OSVs than the 
projected residual (dOSVMP2-P) to recover the same accuracy in the correlation energy. For exam- 
ple, 109 and 105 OSVs respectively yield 99.99% accuracy in the correlation energy for dOSVMP2 
and dOSVMP2-P for the largest peptide [gly]i4; for the (R~20)i9 cluster we require 78 and 73 OSVs 
to reach the same accuracy. 

B. Reaction energies 

Relative energies are the central quantity in chemistry rather than absolute energies. We have 
investigated their accuracy by computing relative energies for some isomerization reactions using 
the different orbital specific virtual approximations. These reactions were selected in Ref. |46| for 
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the good agreement between the canonical MP2 isomerization energies with triple-^ basis sets (cc- 
pVTZ) and the experimental isomerization energies. The results of dOSVMP2-P and OSVMP2 



computations as well as their deviations from the canonical MP2 values are given in Table III The 
dOSVMP2-P computation using only 10% of the OSVs gives a MAD (mean absolute deviation) 
of 2.33 kcal/mol. Using 20% of the OSVs drops the MAD to 0.99 kcal/mol which is chemical 
accuracy. The MAD is further reduced to 0.07 kcal/mol (almost two orders of magnitude) when 
60% of the OSVs are used. 

The MADs of isomerization energies are plotted against the numbers of OSVs in Fig. [3] for 
OSVMP2, dOSVMP2-P and dOSVMP2 schemes. Both dOSVMP2-P and dOSVMP2 display sim- 
ilar accuracies. However, the complete OSVMP2 gives errors that are substantially smaller and 
additionally, these errors decay more rapidly and more smoothly than those of dOSVMP2 and 
dOSVMP2-P, as the number of OSVs used is increased. Nonetheless all methods show a rapid 
decrease in error as the size of the OSV space is increased. 



C. Basis set dependence 



We have investigated the orbital specific virtual orbital dependence of different basis sets using 
the dOSVMP2-P and OSVMP2 approaches. For this we have chosen to use a single glycine molecule 
so that computations with very large basis sets are affordable. As can be seen from Fig. |4j the size 
of the OSV space (required for a given accuracy in the correlation energy) increases much more 
slowly than the size of the underlying basis. Moving from cc-pVDZ to cc-pV5Z |47| . the size of 
the required OSV space for dOSVMP2-P and OSVMP2 increases by a factor of 4-5 while the size 
of the canonical virtual space increases almost by a factor of 10. In fact the size of the OSV space 
needed for a given accuracy appears to increase sub-linearly with the size of the underlying basis. 



D. Visualizing the orbital specific virtual orbitals 

We have visualized the Boys-localized occupied orbitals [IB] and a few associated OSVs for a 



single glycine molecule in Table [TVJ These local occupied HOMO, HOMO-4 and HOMO-8 orbitals 
are chosen to be, respectively, around the N-lone-pair electrons, C-N and C-C bonds along the 
skeleton of glycine. Along each column of Table [TV] each individual OSV exhibits a different shape 



for different occupied orbitals. This can be essentially understood from the definition of the OSVs 



(cf. Eq. (17)) since each OSV has to be adjusted to a particular occupied orbital. For example, 



15 



the LUMO+4 associated to HOMO, HOMO-4 and HOMO-8 demonstrates, respectively, orbital 
locality around the N-lone-pair, C-N and C-C bonds, where the associating occupied orbitals are 
found. 

V. COMPARISON WITH THE PULAY-SAEB0 LOCAL MP2 THEORY 

The Pulay-Saeb0 local correlation approach based on projected atomic orbitals (PAO) is a 
standard against which to compare new approaches to local correlation, such as the orbital specific 
virtual approximations used here. Here we assess both the accuracy and times of the OSVMP2 
and dOSVMP2 approximations relative to the Werner-Schutz formulation of the Pulay-Saeb0 local 
MP2 as implemented in Molpro [6"ll36]. 

A. Potential energy surfaces 

Pulay-Saeb0 local correlation relies on spatial truncation of virtual orbital domains. Discon- 
tinuities on potential energy surfaces (PES) can then arise since the virtual orbital domain size 
defined by spatial truncation is not uniform as the geometry is varied. One prototypical example 
is the propadienone (CH2CCO) molecule that has been recently investigated by Russ et al. [39] 
using local CCSD and MP2. Multiple discontinuities occur in the stretching of the central C=C 
bond of propadienone, even in the vicinity of the equilibrium geometry. Several attempts have 
been made to recover smooth PES in Pulay-Saeb0 theory. By tailoring and fixing virtual do- 
mains [H] the Pulay-Saeb0 local approach can avoid these discontinuities. Explicitly correlated 
R12/F12 methods [SD] reduce the magnitude of discontinuities through the auxiliary excitation 
space [5H [52] . There have also been efforts to use bump functions [53] to smooth discontinuous 
amplitudes. We now reinvestigate this issue using the orbital specific virtual approximations. We 
believe the approach presented here provides a more basic solution. 

The correlation energy PES using different numbers of virtual orbitals (N v ) are presented in 
Fig.[5]for Pulay-Saeb0 PAO MP2, dOSVMP2, dOSVMP2-P and OSVMP2. All computations were 
carried out using the cc-pVDZ basis set jU] with 52 canonical virtual orbitals. For PAO local MP2, 
N v denotes the average size of the pair virtual domain. It can be seen that the PES of the PAO 
based local MP2 with N v = 34 orbitals exhibits three major energy discontinuities in the regions 
of both short and long C=C bonds as well as around the equilibrium C=C bond. When using 
N v = 39 PAOs, five smaller discontinuities in the PAO local MP2 theory appear, ranging from 0.2 
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to 0.4 mEh at 1.269, 1.540, 1.752, 1.978 and 2.352 A. A tiny zig-zag structure can still be seen in 
the vicinity of the equilibrium C=C bond even when using N v = 52 PAOs (see inset of Fig. [5^). 

As for the PES using the orbital specific virtual approximations, both the dOSVMP2, 
dOSVMP2-P and OSVMP2 based curves are smooth even when using much smaller spaces of 
virtual orbitals. The PES of OSVMP2 displays no discernible discontinuities when using 17 OSVs 
(cf. Fig. [Hfj), and when using 22 OSVs the OSVMP2 based PES is already close to the curve 
of canonical MP2. In the case of dOSVMP2 and dOSVMP2-P with N v = 28 OSVs we see very 
tiny breaks (only 0.04 and 0.02 mEh) respectively, at 1.559 and 1.597 A. These discontinuities are, 
nevertheless, one order of magnitude smaller than those of the PAO N v = 34 result discontinuities. 



B. Virtual space size and timings 

The efficiency of both the Pulay-Saeb0 PAO and orbital specific virtual approximations depends 
on the size of the virtual space needed to obtain good agreement with the canonical result. Fig. [6] 
gives the comparison of virtual space sizes between PAO local MP2, dOSVMP2-P and OSVMP2 
needed to recover 99.99% of the correlation energy in [gly] n . It is evident that PAO based local MP2 
needs substantially more virtual orbitals (two times) than dOSVMP2-P and at least four times 
more than OSVMP2. Furthermore, the sizes of the orbital specific virtual spaces for dOSVMP2-P 
and OSVMP2 saturate much more rapidly than those of PAO based MP2 when the molecular size 
increases. The relative advantage of the orbital specific scheme increases as we move to the large 
basis i.e. from cc-pVDZ to cc-pVTZ [47J. 

However, the O 3 scaling term in the orbital specific virtual approximations has a higher prefactor 
than that of the similar term in the Pulay-Saeb0 PAO theory. This affects adversely the CPU times 
to solve the residual equations when O becomes large, even though the orbital specific virtual space 
is much smaller than the PAO space used in the Pulay-Saeb0 scheme. For instance, going from 
[gh/k to [gly]i2, for an accuracy of 99.99% the CPU times of the projective dOSVMP2-P (T s = of 
Table [V]) start to become less favourable than that of the Pulay-Saeb0 theory. The CPU time for 
the full OSVMP2 computations on [gly]i2 is longer than that for the Pulay-Saeb0 implementation 
by a factor of three. 

The computational efficiency of the orbital specific virtual approximations can be greatly im- 
proved, however, by using the screening scheme discussed earlier. We have reinvestigated the CPU 
time to solve the residual equations using screening and the results are presented in Table W\ for 



gly] n and Table VI for polyenes, respectively. The threshold Ts is chosen such that there is only 
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a very minor error in the correlation energy, e.g., the screening errors are restricted to less than 
0.01% of the correlation energy or only a few tenths of kcal/mol in the present study. 

For the longest [gly]i 2 , dOSVMP2-P with T s = 0.20 and OSVMP2 with T s = 0.07 are respec- 
tively speeded up by a factor of 2-3 and 5 compared to those with Ts = 0.00. As a result, both 
screened computations for [gly] n (with a screening error of AE scrn = 0.0001%) are almost two times 
faster than the Pulay-Saeb0 PAO implementation. The screened cc-pVTZ results of dOSVMP2-P 
and OSVMP2 are shown for [gly]4 and [gly]6 in Table |Vj Clearly larger basis sets increase the 
efficiency of screened OSVMP2 computations relative to PAO and dOSVMP2-P. Polyene chains 
are more difficult cases as the correlating orbitals are more extended along the chain than in [gly] n . 
However, if we use a looser accuracy of AE corr < 0.01%, we find that the screened dOSVMP2-P 
and OSVMP2 computations are still faster than the Pulay-Saeb0 implementation by a factor of 2 
and 5, respectively. 

VI. CONCLUSIONS 

In this work we have described the direct orbital specific and full orbital specifc virtual approx- 
imations to local second order M0ller Plesset perturbation theory. These representations of the 
amplitudes have been expressed in a general language of tensor factorization that also encompasses 
many other representations used in electronic structure theory. As we have showed, the orbital 
specific virtual approximation can lead to significant advantages, both in more formal behaviour, 
such as smoothness of potential energy curves, as well as in practical times and accuracies, as com- 
pared to efficient implementations of the local Pulay-Saeb0 correlation ansatz. As for the direct 
versus full orbital specific virtual approximations, when screening is used, the full orbital specific 
virtual approximation is superior. 

There is much to be done along the directions of this work. For example our algorithms are not 
linearly scaling, because we have not investigated efficient representations of the occupied space. 
Furthermore, we expect that significant advantages can be had when applying orbital specific 
virtual type approximations to high body excitations. We conclude by recognizing that the space 
of tensor factorizations is very large, with much remaining to be explored. 
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(b) Matrix Product Decomposition 
FIG. 1: Pictorial representations of (a) Cholesky Decomposition and (b) Matrix Product Factorization. 
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FIG. 2: Pictorial representations of (a) t 



E^*i7*^ (Pulay-Saeb0 PAO approximation), (b) 



tfj = ^Jijitiptijb ( Pair Natural Orbital approximation) and (c) t$ = E^Otf"^ (Orbital Specific 
Virtual approximation) . 
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FIG. 3: MADs (Mean Absolute Deviations) of the isomerization energies as a function of the number of 
OSVs using dOSVMP2, c!OSVMP2-P and OSVMP2. 
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FIG. 4: The number of orbital specific and canonical virtual orbitals using different basis sets (cc-pVXZ, 
X=D, T, Q and 5) for a single glycine molecule. Comparison is made between canonical MP2 (+), 
dOSVMP2-P (unfilled) and OSVMP2 (filled) for accuracies of 99.90% (circle) and 99.99% (square) of the 
correlation energy. Note that OSVMP2 contains the exchange excitations while the dOSVMP2-P does not. 
Consequently, the number of virtual orbitals needed for OSVMP2 is significantly smaller than dOSVMP2, 
roughly half. 
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(a) Pulay-Saeb0's PAO local MP2 




(b) OSVMP2, dOSVMP2 and dOSVMP2-P 

FIG. 5: Potential energy surfaces resulting from (a) Pulay-Saeb0's PAO local MP2 (using Molpro 36J) 
and (b) OSVMP2, dOSVMP2 and dOSVMP2-P for the central C=C bond of propadienone using a cc- 
pVDZ basis set [47]. The equilibrium geometry was obtained by optimizing all internal coordinates of 
propadienone at the MP2/cc-pVDZ level. The displacement of the central C=C bond was 0.001 A and 
other internal coordinates were frozen. The insets magnify the details at the vicinity of the equilibrium 
C=C bonds. 
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FIG. 6: Histograms of virtual orbital numbers for Pulay-Saeb0's PAO local MP2, dOSVMP2-P and OS- 
VMP2 schemes needed to obtain 99.99% accuracy. In the case of PAO local MP2, this number is ob- 
tained from the average pair domain size. Note that OSVMP2 contains the exchange excitations while the 
dOSVMP2-P does not. Consequently, the number of virtual orbitals needed for OSVMP2 is significantly 
smaller than dOSVMP2, roughly half. 
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TABLE I: Comparisons of the number of OSVs needed to obtain different accuracies (99.5%, 99.9% and 
99.99%) in the MP2 correlation energies using dOSVMP2, dOSVMP2-P and OSVMP2 approximations, 
respectively. The cc-pVDZ basis sets [47] were used. The canonical MP2 reference correlation energies were 
obtained using the Molpro program package 36 . N v is the total number of canonical virtual orbitals. 
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TABLE II: Comparison of the relative errors of dOSVMP2, dOSVMP2-P and OSVMP2 correlation energies 
for polyenes using different numbers of virtual orbitals (40, 80, 100 and 140). cc-pVTZ basis sets were used. 
The canonical MP2 reference correlation energies were obtained using the Molpro program package [33] . N v 



is the total number of canonical virtual orbitals. 
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TABLE III: Calculated canonical MP2, dOSVMP2-P and OSVMP2 isomerization reaction energies 
(kcal/mol). The reactions and corresponding molecular geometries were taken from [46] The dOSVMP2-P 
and OSVMP2 results are given as the first and second row for each reaction. The errors given in the paren- 
theses are the deviations relative to canonical MP2 reaction energies. The Ahlrichs-TZV [Ml [55] basis sets 
augmented by the cc-pVTZ polarization functions (2dlf, 2pld) [47] were used together with based on the 
frozen-core approximation. The percentage is the first row indicating the fraction (%) of the orbital specific 



virtual space used in the calculation. 



J. Vj O C*i(_. Li 1 U 1 1 o 




MP 2 




10% 






20% 






40% 






60% 




7 




9.26 


10 


75 


1 


49) 


11 


1 


' I 


84) 


10 


06 


f 


80) 


9 


46 




20) 








10 


54 


( 1 


28) 


9 


76 


( () 


50) 


9 


28 


( o 


02) 


9 


26 


( o 


00) 


8 




22.20 


20 


63 


'_l 


57) 


21 


30 


(-[) 


90) 


21 


82 


f-D 


38) 


22 


1 


u 


10) 








21 


29 


[-0 


91) 


22 


00 


(-0 


20) 


22 


19 


(-0 


01) 


22 


20 


( o 


00) 


9 




6.96 


5 


17 


(-1 


79) 


5 


74 


(-1 


22) 


6 


92 


(-0 


04) 


6 


94 


(-0 


02) 








C 




54 


(-0 


42) 


C 




91 


(-0 


05) 


6 


96 


( o 


00) 





96 


( o 


00) 


1Z 




47.13 


6 1 


61 


(-9 


OZ ) 


4 ( 


25 


( o 


12) 


4 i 


01 


(-0 


1 o\ 
12) 


4 1 


09 


(-0 


U4J 








4 r i 

-to 


76 


(-1 


SI) 


4-7 


03 


(-0 




47 


12 


(-0 




47 


12 


( o 


UU 1 


la 




11.52 


1 n 

1U 


75 


(-0 


77) 


11 


13 


(-0 


3y ) 


11 


53 


( 


01) 


11 


52 


( o 


UU ) 








1 1 


28 


(-0 


Z1) 




48 


(-0 




1 1 

_L L 


52 


( o 


UU ( 


1 1 


52 


( o 


UU J 


zl 




1.08 


1 


27 


( o 


iy ) 


1 


03 


(-0 


(Jo ) 


1 


11 


( o 


U3 ) 


1 


07 


(-0 


Ul) 








1 


15 


( o 


07) 


1 


09 


( o 


01) 


1 


08 


( o 


00) 


1 


08 


( o 


00) 


24 




12.56 


13 


44 


( o 


88) 


12 


85 


( o 


29) 


12 


36 


(-0 


20) 


12 


48 


(-0 


08) 








12 


68 


( o 


12) 


12 


52 


(-0 


04) 


12 


56 


( o 


00) 


12 


56 


( o 


00) 


28 




31.12 


33 


90 


( 2 


78) 


33 


10 


( 1 


98) 


31 


74 


( o 


62) 


31 


32 


( 


20) 








32 


64 


( 1 


52) 


31 


45 


( o 


33) 


31 


14 


( o 


02) 


31 


12 


( o 


00) 


32 




7.35 


3 


04 


(-4 


31) 


5 


70 


(-1 


65) 


7 


01 


(-0 


34) 


7 


29 


(-0 


06) 








6 


20 


(-1 


15) 


7 


15 


(-0 


20) 


7 


35 


( o 


00) 


7 


35 


( o 


00) 


34 




6.98 


7 


00 


( o 


02) 


5 


52 


(-1 


46) 


6 


71 


(-0 


27) 


6 


97 


(-0 


01) 








6 


62 


(-0 


36) 


6 


96 


(-0 


02) 


6 


98 


( o 


00) 


6 


98 


( o 


00) 


Mean Absolute Deviation 


(kcal/mol) 




2. 


13 






0. 


99 






0. 


28 






0. 


07 












0. 


74 






0. 


15 






0. 


01 






0. 


JO 




Mean Deviation 


(kcal/mol) 




-1 


26 






-0 


14 






0. 


Jl 






0. 


01 












-0 


15 






0. 


02 






0. 


JO 






0. 


00 




Maximum Deviation 


(kcal/mol) 




9. 


52 






1. 


98 






0. 


80 






0. 


20 












1. 


52 






0. 


50 






0. 


J2 






0. 


00 





30 



TABLE IV: Contour plots of a few Boys-localized occupied orbitals (HOMO, HOMO-4 and HOMO-8) and 
associated OSVs (LUMO, LUMO+4 and LUMO+5) for an isolated single glycine molecule. The OSVs are 



sorted in descending order according to the singular values (cf. Eq. (18)). HOMO and LUMO correspond 
to the orbital, respectively, that has the highest occupied orbital energy and the largest singular value. 
HOMO-4 and HOMO-8 give the fourth and eighth localized occupied orbitals with orbital energies below 
the HOMO. LUMO+4 and LUMO+5 are the fourth and fifth OSVs with the singular values above the 
LUMO. 



k 


Local Occ. 


LUMO 


LUMO+4 


LUMO+5 


HOMO 






k 




HOMO-4 






k 




HOMO-8 
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TABLE V: CPU times (i solv , seconds) to solve the dOSVMP2-P and OSVMP2 residual equations, using 
different screening thresholds Ts for [gly]„ chains with n = 4,6,8,10,12. 103 and 46 OSVs have been 
used respectively for dOSVMP2-P and OSVMP2 in order to obtain accuracies of 99.99% for all T s = 0.00 
computations. AE scrn is the additional relative error in the correlation energy introduced by screening 
Tg > 0.00 as compared to canonical MP2. The cc-pVDZ basis set [?7] (denoted as D) was used for all 
computations and results with the cc-pVTZ basis set (denoted as T) are also reported for [gly]4 and [gly]6 
molecules. With the cc-pVTZ basis set N v = 235 ([gly] 4 ) and N v = 243 ([gly] 6 ) were used for dOSVMP2-P 
while N v = 92 ([gly] 4 ) and N v = 94 ([gly] 6 ) were used for OSVMP2. 







N v = 103 for dOSVMP2-P 


N v = 46 for OSVMP2 




PAO 


T s = 0.00 


T s - 


= 0.20 


T S = 


= 0.40 


T s = 0.00 


T s - 


= 0.07 


T s = 


= 0.20 


[giy]„ 


^solv 


^solv 


^solv 


AE scrn 


^-solv 


AE 

' 1 1 'sera 


^solv 


^solv 


AE 


^solv 


AE 


4(D) 


400.8 


470.1 


378.1 


0.0001% 


340.3 


0.0002% 


624.7 


307.4 


0.0001% 


242.5 


0.008% 


4(T) 


4819.4 


5374.8 


4041.2 


0.0001% 


3614.2 


0.0001% 


4270.1 


1691.6 


0.0005% 


1414.5 


0.01% 


6(D) 


1280.6 


1303.5 


832.5 


0.0001% 


749.2 


0.0004% 


2104.1 


741.7 


0.0001% 


568.9 


0.009% 


6(T) 


12767 


16429 


10032 


0.0002% 


8872.6 


0.0004% 


13911 


4015.3 


0.0005% 


3328.2 


0.01% 


8(D) 


2436.7 


2748.8 


1465.7 


0.0001% 


1309.7 


0.0004% 


4739.4 


1303.8 


0.0001% 


984.8 


0.010% 


10 (D) 


4038.1 


5050.7 


2283.7 


0.0001% 


2027.1 


0.0005% 


8981.6 


2023.5 


0.0001% 


1515.1 


0.011% 


12 (D) 


5147.0 


8248.9 


3270.5 


0.0001% 


2912.2 


0.0005% 


15219.3 


2901.4 


0.0001% 


2153.8 


0.011% 
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TABLE VI: CPU times (t solv , seconds) to solve the dOSVMP2-P and OSVMP2 residual equations, using 
different screening thresholds Tg, for polyenes. The PAO computation recovers about 99.26% of the cor- 
relation energy. 100 and 40 OSVs have been used for dOSVMP2-P and OSVMP2, respectively AE scrn is 
the additional relative error of correlation energy introduced by screening T s > 0.00 compared to canonical 
MP2. For Ts = 0.00 the percentage (AE corr ) of the canonical MP2 energy recovered is also reported. The 
cc-pVTZ basis set was used for all computations. 



N v = 100 for dOSVMP2-P N v = 40 for OSVMP2 





PAO 


Ts 


= 0.00 
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T s 


= 0.20 
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^•solv 


^solv 


AE 
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AE 
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AE 
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AE 
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AE 




20.4 


26.2 


99.61% 


24.3 


0.001% 


23.1 


0.004% 


15.2 


99.54% 


11.2 


0.002% 


CsHio 


53.5 


50.4 


99.48% 


43.3 


0.002% 


40.2 


0.007% 


33.0 


99.46% 


19.6 


0.004% 


C10H12 


110.3 


86.2 


99.39% 


68.8 


0.003% 


63.4 


0.009% 


62.8 


99.41% 


31.8 


0.005% 


C12H14 


211.5 


135.0 


99.32% 


98.6 


0.004% 


89.6 


0.009% 


103.4 


99.37% 


44.3 


0.007% 


C14H16 


313.8 


224.4 


99.27% 


137.1 


0.004% 


123.2 


0.009% 


163.7 


99.34% 


62.1 


0.008% 
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